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Abstract 



In air shower experiments information about the initial cosmic ray particle or about the shower 
development is obtained by exploiting the correlations between the quantities of interest and the 
directly measurable quantities. It is shown how these correlations are properly treated in order to 
Q> ■ obtain unbiased results. As an example, the measurement of the average penetration depth as a 
, function of the shower energy is presented. 



t^)! 1 Introduction 

In air shower experiments cosmic ray particles are studied in an indirect way: the distributions of 



< 



the interesting quantities (Xorig) of the initial cosmic ray particle (like its nature and energy) or of 
the air shower (like the penetration depth) have to be inferred from the distributions of measurable 
quantities (Xmcas) (like particle and C-light densities at detector level). Using Monte Carlo (MC) 
^ I simulations, in which the interaction of the cosmic ray particle with the atmosphere, the shower devel- 
opment and also the properties of the detector are simulated, one is able to establish the correlations 
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Q , between the measurable quantities Xmeas and the interesting quantities Xorig. These correlations are 
i then used to determine the distributions of Xorig from the experimental distributions of Xmeas- 
I The aim of this paper is to demonstrate that it is essential to treat the correlations in a mathe- 

■ matically correct way in order to avoid biases in the results. As an example the determination of the 
average penetration depth Xmax of air showers in the atmosphere as a function of the shower energy 
Qh| E is presented, using data from the HEGRA array of scintillator and wide-angle C-light detectors. 
Q ■ In this example, the number of shower particles (Ng) and the C-light radius at detector level 

(Rl) are taken as measurable quantities. Ng is determined from the particle densities as measured 
^ I by the matrix of scintillation detectors and Rl is obtained as the inverse of the slope of the lateral 
C-light distribution as measured by the matrix of C-light detectors. The quantities of interest ("true" 
quantities) are the shower energy E and the penetration depth Xmax* 

In first approximation, Ng carries mainly information on E, and Rl mainly on Xmax- While the 
I correlation between Rl and Xmax is quite independent of the nature (or atomic number A) of the 
cosmic ray particle (Lindner, 1998), this is not the case for the correlation between Ng and E. It has 
been shown in (Cortina, 1997) that by using a modified Ng (N^ = Ng • Rl°, with a depending on the 
zenith angle) the correlation between Ng and E is quite independent of A. 

2 The Method 

The task is now to determine the 2-dimensional distribution of the variables (log E, Xmax) from 
the 2-dimensional distribution of the variables (log Ng, 1/Rl). The simplest way of doing this is to 
use the average correlations 

< log >= fi(log E); < 1/Rl >= f2(Xmax) 

as determined from a sample of Monte Carlo events. In this procedure a one-to-one correlation is 
assumed between log E and log and between Xmax and 1/Rl respectively. In addition, a possible 
Xmax dependence of the log E — log Ng correlation, and a log E dependence of the Xmax — 1/Rl 



correlation is ignored. This procedure will later be referred to as the "one-to-one correlation proce- 
dure" . 

A mathematically correct approach is to use the full correlation between (log E, Xmax) and 
(log Ng, 1/Rl)- It is convenient to define a grid in the (log E, Xmax) plane and a grid in the (log Ng, 
1/Rl) plane, and some numberings (i = 1 to Nxy) and (j = 1 to Nuv) of the resulting bins in the 
two planes respectively. The bin content of bin (j, i) of the distribution (log Ng, 1/Rl) versus (log E, 
Xmax) for the sample of Monte Carlo events may be denoted by Gji. The full correlation can then be 
written as 

«' = erfe 

gji^ describes how a particular pair of values (log E, Xmax), defined by a specific value io of i, is 
transformed into a distribution of (log N^jI/Rl), given by gji^ (j = 1 to Nuv)- The division by 
{J2k Gki) was done to make gji independent of the (log E, X^ax) distribution in the Monte Carlo 
sample. 

If the experimental distribution of (log Ng, 1/Rl) is denoted by aj (j = 1 to Nuv) and the distribu- 
tion to be determined in (log E, X^ax) by bi (i = 1 to Nxy), then bi has to fulfill the condition 

aj=E(gji-bi) (2) 

The condition ^ ensures that the full correlations between the measured and the true quantities are 
taken into account. 

Determining the distribution bi from a measured distribution aj, with known response matrix gji, 
is a typical unfolding problem. The main point of the unfolding methods is to impose, in addition to 
((2), certain smoothness conditions on the distribution bi in order to avoid strong fluctuations of bi, 
which arise from statistical fluctuations of aj . In the example discussed here the method of reduced 
crossed entropy (MRX) is applied (Schmelling, 1994). In this method, a kind of smoothness condition 
is imposed by requiring the solution bi to be close to a prior distribution b?"^'™. b?"^'™ may be some 
guess of the true distribution. Usually the result bi is quite independent of the choice of b?"^'™, so 
that b?"^'™ may be set to a constant. 

3 Results 

The response matrix gji for the example discussed here is shown in Fig. 1. The ordinate corre- 
sponds to bins in the (log E, X^ax) plane, the abscissa to bins in the (log Ng, 1/Rl) plane, gji was 
obtained by averaging the response matrices for different chemical elements (A = 1, 4, 16 and 56) 
and by smoothing the average response matrix in the following way: gji was parametrized as a 
2-dimensional Gaussian distribution in the variables (u, v) = (log Ng, 1/Rl) 
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where the 5 parameters p = u, v, cju, cr^ and p where assumed to be linear functions of log E and Xj, 
(3 parameters for each of the 5 parameters p). The 5 x 3 = 15 free parameters were determined by 
fitting expression Q to the average response matrix. The fitted values of the parameters characterize 
in detail the behaviour of the correlations between (log E, X^ax) and (log Ng,l/RL) and their de- 
pendence on log E and Xmax- In particular one finds: In very good approximation, < 1/Rl > (= v) 
is only a function of Xmax- < log Ng > (= u) is mainly a function of log E, with some additional 
dependence on Xmax- The parameter p, which describes the correlation between log Ng and 1/Rl at 
fixed (log E, Xmax), is a function of log E: at low log E log Ng and 1/Rl are anti-correlated, whereas 
they are positively correlated at higher log E. All these properties of the correlations are, of course, 
taken into account in the unfolding procedure. 



The experimental distribution of (log Ng,l/RL) is shown in Fig. 2 (Kornmayer, 1999). It should 
be noted that the measurements presented in this figure are based on preliminary data and are not 
the final official HEGRA results. By applying the MRX method a distribution of (log E, Xmax) is 
obtained ("unfolded" distribution) which is displayed in Fig. 3. Forming the average Xmax for each 
bin of log E yields the result for the elongation plot < Xmax > versus log E, shown in Fig. 4a. In 
Fig. 4b the RMS of X^ax is plotted as a function of log E. 

For comparison, in Fig. 4 the results from the "one-to-one correlation procedure" (see above) 
are also plotted. It can be seen that the latter procedure underestimates < X^ax > by ~ 30 g/cm^ 
at low log E and by ~ 10 g/cm^ at high log E. The points in the bin of highest energy should be 
taken with care because they are based on low statistics both in the experimental data and in the 
MC sample. No systematic differences between the two methods are seen for the RMS of X^ax (Fig- 
4b). Knowing that the one-to-one correlation procedure yields biased results one can try to correct 
the results by applying additional correction factors to < X^ax >; which are determined from MC 
events. However, these correction factors will in general depend on the details of the MC simulation, 
in particular on the distribution of (log E, Xmax)- 

By construction, the result of the unfolding procedure does not depend on the underlying MC 
distribution of (log E, Xmax)- By fulfilling the condition (jSJ (at least approximately, see Fig. 2), 
it takes into account the full correlations between the measured and the true quantities. Of course, 
these correlations and thus also the result for the elongation plot will depend on the model used in 
the MC simulation. How they depend on the MC model can be studied by doing the unfolding for 
different response matrices, corresponding to different MC models. 

A study of the dependence of the results from the one-to-one correlation procedure on the MC 
model will be less conclusive because effects due to differences between the MC models and effects 
due to using a mathematical incorrect procedure are not well separated. 

Since the average penetration depth of air showers depends on the nuclear mass number A of the 
cosmic ray particle inducing the air shower, a measurement of < Xmax > as a function of E can be 
used to obtain information on the chemical composition of cosmic rays (see for example Roehring, 
1999). 

It should however be noted that from the same experimental data information about the chemical 
composition can also be obtained in a more direct way: one possibility is to start from the experi- 
mental 2-dimensional distribution of (log Ng, 1/Rl) and apply the unfolding procedure to obtain the 
2-dimensional distribution of (log E, log A). In this case the response matrix would explicitly depend 
on log E and log A and one would not have to rely on an A-independence of the response matrix, 
as was the case for the example discussed in this paper. The A-independence was necessary in order 
to determine the distribution of the penetration depth. If one is only interested in the chemical 
composition a knowledge of this distribution is not required and the 2-dimensional distribution of 
(log E, log A), which contains all the information about the chemical composition as a function of E, 
can be obtained directly. 
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Figure 1: Response matrix gji. The abscissa cor- 
responds to bins in the (log Ng, 1/Rl) plane: In 8 
consecutive bins 1/Rl increases from 0.0 to 0.035 
m~^. Each block of 8 consecutive bins is for one 
bin in log Ng. In 8 consecutive blocks log Ng in- 
creases from 7.27 to 10.07. The ordinate corre- 
sponds to bins in the (log E, Xmax) plane: In 8 
consecutive bins Xmax increases from 320 to 800 
g/cm^. Each block of 8 consecutive bins is for 
one bin in log E. In 8 consecutive blocks log E 
increases from 1.5 to 4.3. 
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Figure 2: Experimental distribution of 1/Rl in 
different bins of log Ng (full circles). The his- 
togram represents the result from applying the 
response matrix to the unfolded distribution of 
(log E, Xmax)- The good agreement between the 
two distributions shows that relation (2) between 
the measured and unfolded distribution is well ful- 
filled. The abscissa in this figure corresponds to 
the abscissa in Fig. ^ 
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Figure 3: Unfolded distribution of X^ax in differ- 
ent bins of log E (open circles). The histogram 
represents the prior distribution used in the un- 
folding procedure. The abscissa in this figure cor- 
responds to the ordinate in Fig. ^ 
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Figure 4: a) Average penetration depth < X^ax > 
and b) RMS of X 

max as functions of E. The full 
circles represent the results of the unfolding pro- 
cedure using the full response matrix. If one- 
to-one correlations are assumed between 1/Rl 
and Xmax) and log Ng and log E respectively, the 
points represented by open circles are obtained 
("diagonal response matrix"). While both meth- 
ods yield consistent results for RMS of Xmax, the 
one-to-one correlation procedure in general un- 
derestimates < Xmax > by 10 to 30 g/cm^. 



